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The interlayer sHding energy landscape of hexagonal boron nitride (/i-BN) is investigated via a 
van der Waals corrected density functional theory approach. It is found that the main role of the van 
der Waals forces is to "anchor" the layers at a fixed distance, whereas the electrostatic forces dictate 
the optimal stacking mode and the interlayer sliding energy. A nearly free-sliding path is identified, 
along which bandgap modulations of --^ 0.6 eV are obtained. We propose a simple geometrical model 
that quantifies the registry matching between the layers and captures the essence of the corrugated 
/i-BN interlayer energy landscape. The simplicity of this phenomenological model opens the way to 
the modeling of complex layered structures, such as carbon and boron nitride nanotubes. 
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The interlayer potential landscape in layered materials 
is essential for understanding their mechanical and elec- 
tromechanical behavior. For example, nanoelectrome- 
chanical systems (NEMS) based on low dimensional 
structures of such materials often rely on mechanical 
deformations such as twisting [ll-Q and bending 
These processes involve relative sliding of the layers, 
which exhibits a corrugated energy landscape even in 
atomically flat systems, such as graphite and hexago- 
nal boron nitride (/i-BN). This corrugation arises from 
the non-uniform charge density distribution around the 
atomic positions within each layer 0-01 • It is well ac- 
cepted that the most important factors that govern cor- 
rugation in such system are electrostatic and dispersion 
interactions. However, a clear picture of how the com- 
plex interplay between these factors determines the cor- 
rugated energy landscape and its manifestation in unique 
material properties has not emerged yet. 

Previous efforts towards the understanding of these 
phenomena have utilized density functional theory 
(DFT) with (semi-)local approximations [ol-fl^. How- 
ever, these methods do not provide an appropriate de- 
scription of dispersive interactions. Several approaches 
within DFT have been developed to overcome this prob- 
lem, demonstrating successful treatment of dispersion ef- 
fects in graphite 13-15| and /i-BN [l3, 14, 1^. However, 
to the best of our knowledge, such studies have not ad- 
dressed the question of corrugation. 

In this letter, we present a theoretical study of the 
intricate interplay between dispersion and electrostatic 
interactions in layered materials. As an illustrative ex- 
ample we choose the case of /i-BN, where both types of 
interactions arc expected to have a considerable influ- 
ence on the landscape of the interlayer potential. Using 
a first-principles van der Waals (vdW) corrected DFT ap- 
proach, we show that dispersion interactions play the role 
of fixing the interlayer distance, while the electrostatic 
forces determine the optimal stacking mode and the in- 



terlayer sliding corrugation. In addition, we predict the 
existence of a nearly free-sliding path along which con- 
siderable bandgap modulations are obtained. Finally, we 
propose a simple geometrical model that quantifies the 
registry matching between the layers and captures all the 
important physical features appearing in the corrugation 
energy landscape of /i-BN. 
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FIG. 1: Binding energy curves of bulk and bilayer /t-BN at the 
AA' stacking mode, calculated using the PBE functional with 
and without the TS-vdW correction. Energies are presented 
with respect to two infinitely separated BN sheets. 

We start by demonstrating the importance of appro- 
priate treatment of dispersion interactions when mod- 
eling layered materials in general and /i-BN in particu- 
lar. To this end, we calculate the interlayer distance de- 
pendence of the binding energy (BE) of /i-BN using the 
generalized gradient approximation (GGA) of Perdew, 
Burke, and Ernzerhof (PBE) [l3| augmented with the 
Tkatchenko-Scheffler vdW (TS-vdW) correction 
Within this approach, the leading order pairwise Cq/R^ 
dispersion correction is added to the internuclear energy 
term, where the Cg coefficients and the vdW radii are 
determined directly from the DFT ground state electron 
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density and reference values for the free atoms. 

The geometry of a single BN sheet was optimized 
with the PBE functional using two-dimensional periodic 
boundary conditions, |20| as implemented in the Gaus- 
sian suite of programs [21| using the large quadruple- 
zeta Weigend-Ahlrichs basis set (Def2QZVP) [23|. Unit 
cells of bilayer and bulk /i-BN were then constructed 
by stacking single BN sheets in the experimentally ob- 
served AA' stacking mode (sec Fig. ^jp). The BE of 
these structures was calculated as a function of the in- 
terlayer distance at fixed intralayer coordinates with and 
without the TS-vdW correction, as implemented in the 
FHI-aims code We used the tier-2 numerical atomic- 
centered orbitals basis set, known to yield binding ener- 
gies converged to the meV/atom level [31, verified here 
by comparison with tier-3 basis set calculations for se- 
lected stacking modes. 

As shown in Fig. [TJ the PBE+vdW equilibrium inter- 
layer distance calculated for bulk /i-BN (solid blue line) 
is 3.33 A, in perfect agreement with the experimental 



value [2J]. This is in stark contrast to the result ob- 



tained from the pure PBE functional: an overestimated 
equilibrium inter layer distance of 4.17 A accompanied 
by a ver y in odest BE of 2.05 meV/atom (dashed pur- 
ple line) |25l . [26j . The lattice constant obtained with the 
PBE-|-vdW approach is consistent with the value of 3.31 
A [3l obtained previously using the adiabatic-connection 
fluctuation-dissipation theorem, but here the computa- 
tional cost is significantly reduced. The PBE-|-vdW bind- 
ing energy for bulk /i-BN is 85.9 mcV/atom, which is 
somewhat higher then the BE of 56 mcV/atom recently 
reported for graphite flif . 

Similarly, in the case of bilayer /i-BN, PBE yields a 
small BE of 1.0 meV/atom and an overestimated inter- 
layer distance of 4.22 A, whereas PBE-|-vdW predicts a 
much larger BE of 38.1 meV/atom and an interlayer dis- 
tance of 3.37 A [23|. We note that the interlayer binding 
in bilayer /i-BN is weaker than that obtained for bulk h- 
BN. This is due to the fact that in the bulk system each 
layer interacts with two adjacent layers, instead of one, 
as well as with additional layers that are farther away. 
We note that Rydberg et al. [31, reported a nearly in- 
distinguishable behavior of bulk and bilayer /i-BN. How- 
ever, they obtained an overestimated interlayer distance 
of 3.63 A. 

Having established the validity of PBE-|~vdW approach 
for the description of the interlayer coupling in /i-BN, we 
now apply it to the study of corrugation in the bilayer 
system. Starting from the AA' stacking mode with an 
interlayer distance of 3.37 A, we perform a set of lateral 
shifts of one /i-BN layer parallel to the basal plane of of 
the other. At each shifted configuration we calculate the 
total energy of the bilayer system. The resulting sliding 
energy landscape is presented in Fig. [2^. 

In order to quantify the role of vdW interactions for 
the sliding process, we compare the changes in total cn- 




FIG. 2: Corrugation of bilayer /i-BN as a function of lateral 
interlayer shifts: (a) Energy surface calculated with respect to 
the AA' stacking mode using the PBE-fvdW approach. The 
yellow line indicates the nearly-free-sliding path (see text); 
(b) Sliding energy along the nearly-free-sliding path. Insets 
are illustrations of selected stacking configurations along the 
path; (c) Registry index surface. Inset: Illustration of the 
geometric model and the dilferent overlap terms. 



ergies with the changes in the vdW contribution upon 
interlayer sliding. We find that the maximal total energy 
change obtained is 26 meV/unit-cell while the largest 
change in the vdW correction is 5 meV/unit-cell (not 
shown). Comparing with the vdW contribution to the 
BE, which is two orders of magnitude larger (see Fig.[T|), 
we conclude that the main role of the vdW interactions is 
to "anchor" the layers at the appropriate interlayer dis- 
tance. Once stacking is established, the sliding energy 
profile at fixed interlayer distance is governed by electro- 
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static interactions resulting from the polar nature of the 
B-N bond 0, Q ■ These conclusions are further supported 
by the fact that the general shape of the sliding energy 
surface remains qualitatively unchanged when perform- 
ing the calculation without the vdW correction, provided 
that the layers are kept at the same interlayer distance. 

Apart from the stable AA' configuration, two other 
high symmetry configuration are obtained. We mark 
them as ABi (See Fig. [2}d) and AB2 in accordance with 
the standard nomenclature of graphite stacking. In the 
ABi mode, the boron atoms of the top layer lie exactly 
above those of the bottom layer, while the nitrogen atoms 
are located above the center of a BN hexagon. The AB2 
mode is similar but with interchanged positions of the 
boron and nitrogen atoms. Despite the similar geometry, 
the nature of the two stacking modes is surprisingly dif- 
ferent 0. The ABi mode is stable with a total energy 
comparable to that of the AA' mode, whereas the AB2 
mode is unstable with a total energy 26 meV/unit-cell 
higher than the AA' configuration. We note that with 
the pure PBE functional the AA' stacking mode is the 
global minimum of the sliding energy landscape, lower by 
only ~ 3.5 meV/unit-cell than the ABi mode. When ap- 
plying the vdW correction, the ABi stacking mode turns 
out to be lower by ^ 1.4 mcV/unit-cell than the AA' 
configuration. Because this difference is within the ac- 
curacy limits of our approximations we do not attribute 
any physical importance to it. 

An interesting feature that appears in the sliding en- 
ergy surface is the existence of a nearly free-sliding path, 
indicated by the yellow line appearing in Fig. [2^1. Fig.[2lD 
presents the total energy differences with respect to the 
AA' stacking mode along this path, showing total energy 
variations of ~ 2 mcV/unit cell. This suggests the /i-BN 
will exhibit highly anisotropic mechanical and electrome- 
chanical properties. 

The calculated corrugation energies presented above 
indicate that different stacking modes exhibit consider- 
ably different energetic stability. Due to the inverse cor- 
relation between stability and bandgap, one may expect 
that upon interlayer sliding the electronic properties of 
layered /i-BN will vary. In order to quantify this elec- 
tromechanical effect we calculate the bandgap of bilayer 
/i-BN as a function of interlayer sliding. We perform these 
calculations using the screened-exchange hybrid approx- 
imation of Hcyd, Scuscria, and Ernzcrhof (HSE06) pll . 
[28| . which is expected to reproduce experimental optical 
bandgaps of bulk semi-conductors [29| , and the double- 
C polarized 6-31G** Gaussian basis set [s^l- This ap- 
proach has been recently shown to describe the physical 
properties of graphene-based materials with exceptional 
success [3ll |. 

The sliding induced bandgap variations are presented 
in Fig. [3l A steep 0.6 eV bandgap decrease is obtained 
as the layers are laterally shifted away from the AA' con- 
figuration. This change is about 10% of the calculated 




FIG. 3: Bilayer-/i-BN bandgap variations as a function of 
lateral interlayer shifts at a fixed interlayer distance of 3.37 

A. 

bandgap value of 6.05 eV for the AA' stacking mode, in- 
dicating a strong electromechanical response of bilayer 
/i-BN towards interlayer sliding. Furthermore, we find 
that in addition to the large variations in magnitude, the 
bandgap changes its nature from direct to indirect dur- 
ing the interlayer sliding process, in agreement with the 
calculations of Liu ct al. [9| . It is interesting to note that 
large bandgap variations are obtained along the nearly 
free-sliding path, meaning that the electronic properties 
of /i-BN will be sensitive to stress applied along this di- 
rection. 

From what we have presented thus far, it is clear that 
there is an intimate relation between the registry match- 
ing of the different 2D crystalline sheets and the energetic 
stability of the various stacking modes in layered mate- 
rials. Previous studies have addressed this issue only 
in qualitative terms regarding different configurations as 
having "good" or "bad" stacking 

Here, we present a simple and intuitive model that 
quantifies the registry matching using basic geometric 
considerations. As stated above, the main contribution 
to the corrugation energy in /i-BN comes from the elec- 
trostatic interactions between the atomic sites, which 
carry apartial charge due to the polar nature of the B-N 
bond [71, y. In order to mimic this effect, we ascribe to 
each atom in the unit cell a circle centered around its 
position. Considering the projection on a plane parallel 
to the layers, as illustrated in Fig.[2j:, we mark by Sij the 
overlaps between the circle centered around the i atom 
belonging to the top layer and the j atom of the bottom 
layer {i and j being either N or B). We now define the 
registry index as: 

_ [Snn — '5'ivAr ) + {Sbb — S^g ) — [Snb — S'j^B ) 

~ / qAA _ qAA' \ \ I oAA _ qAA' \ _ ( qAA qAA' \ ' 
WWJV ^NNJ^y^BB '^BB ) N B N B I 

where the factors Sfj^ and Sfj"^ are the respective over- 
laps at the AA' and AA stacking modes, introduced for 
normalization purposes. Here, AA denotes the case were 
the two layers are completely eclipsed and in perfect reg- 
istry. With this definition, the registry index, i?, is lim- 
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ited to the interval [0, 1] where the minimal (maximal) 
value is obtained for the AA' (AA) mode. 

It is now possible to plot the registry index as a func- 
tion of the lateral interlayer shifts while using the ratio 
between the radii of the N (r^v) and B (rs) circles as 
a single fitting parameter to obtain good agreement be- 
tween the registry index and the corrugation energy sur- 
faces. In Fig. [2jl we plot the registry index surface for 
rjv = 0.50Rbn and = O.WRbn, where Rbn ~ 1-45 A 
is the equilibrium BN bond length in /i-BN. This ratio be- 
tween the two radii takes into account the non-uniform 
charge distributions around the B and N atomic sites 
where the nitrogen (boron) atom has a larger (smaller) 
electron cloud around it. Clearly, a good agreement is ob- 
tained between the registry index calculated via the sim- 
ple geometric model and the corrugation energy surface 
obtained from state-of-the art DFT calculations. There- 
fore, we conclude that it is possible to capture all the 
important physical parameters that govern the corruga- 
tion process using simple geometric considerations [32|. 

In summary, we have studied the corrugated sliding en- 
ergy landscape in layered /i-BN via a van der Waals cor- 
rected DFT approach. The delicate interplay between 
different contributions was investigated, revealing that 
the main role of the van der Waals forces is to anchor 
the layers at a fixed distance, whereas the electrostatic 
forces dictate the optimal stacking mode and the in- 
terlayer sliding corrugation. A nearly free-sliding path 
has been identified, along which bandgap modulations of 
^ 0.6 eV are predicted to occur. A simple geometrical 
model that quantifies the registry matching between the 
layers and captures the essence of /i-BN interlayer sliding 
was presented. The model is not limited to the case of 
/i-BN and can be easily extended to treat other layered 
materials, such as graphite, and generalized to describe 
more complex structures such as niultilaycrcd nanotubes 
of different types. Following the guidelines presented in 
this letter, it is possible to identify stable configurations 
and quantify the corrugation in complex layered systems 
based on simple geometrical arguments with a negligible 
computational effort. 
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